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A numerical method of calculating the non-Markovian evolution of a driven atom radiating into 
a structured continuum is developed. The formal solution for the atomic reduced density matrix is 
written as a Markovian algorithm by introducing a set of additional, virtual density matrices which 
follow, to the level of approximation of the algorithm, all the possible trajectories of the photons 
in the electromagnetic field. The technique is perturbative in the sense that more virtual density 
matrices are required as the product of the effective memory time and the effective coupling strength 
become larger. The number of density matrices required is given by 3 M where M is the number 
of timesteps per memory time. The technique is applied to the problem of a driven two-level atom 
radiating close to a photonic band gap and the steady-state correlation function of the atom is 
calculated. 
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I. INTRODUCTION 



Photonic band-gap (PBG) materials give rise to an altered density of states (DOS) for the electromagnetic field, 
which leads to a fundamentally different atom-radiation field interaction [Q-|8|. When atoms emit light into these 
structures, they exhibit inhibition of spontaneous emission |^,[n^|, localization of light, and the formation of an atom- 
photon bound state ]TT|-|l3|]. These effects are due to entanglement between the atomic and photonic states, and as 
such represent a failure of the Born-Markov approximation, which allows the derivation of a master equation for the 
reduced density matrix of the atom and has been an extraordinarily useful tool in the field of theoretical quantum 
optics for treating spontaneous emission. The spontaneous emission of photons into a PBG has been described 
analytically, but describing such a system with the added noniinearity of driving has remained a theoretical challenge 
HQ . This paper outlines a numerical method that calculates the resonance fluorescence of a driven atom in a PBG 
Oh, material, ft may be applicable to other weakly non-Markovian systems. 

+Jj ■ A photonic band gap is a range of frequencies in which no photon propagation can occur. An atom near one of 
these band edges will experience an interaction with an altered electromagnetic vacuum, and will therefore decay in 
a different manner. Trivially, there will be no emission of photons into the band gap, but the altered DOS near the 
band edge will have further effects on the behavior. The case of spontaneous emission was treated by John and Quang 
JTo| who, due to the restricted state space in this case, were able to solve the equations of motion for the system 
analytically by Laplace transform methods. They found that the atom would not decay to the ground state due to 
the presence of a localized atom-photon bound state which was more strongly bound with stronger coupling. Near 
a band edge, the energy of the propagating photons varies quadratically with the wavevector in the same way as a 
massive particle propagating in free space. The model of this system is therefore virtually identical to that of a single 
mode atom laser |lq-|lq1 . 

The extension of the work of John and Quang to a driven atom has remained a theoretical challenge. When 
pumping was added to the atom laser model, it allowed access to a potentially infinite ladder of states within the 
lasing mode. Fortunately, the linearity of the problem allowed a closed set of equations for the relevant expectation 
values to be developed, and the steady state output spectrum of the system was calculated [^9|. In contrast, a driven 
two level atom in a PBG material has a much smaller system (two states), but the nonlinear behavior of the atom 
due to the driving, in conjuction with the output coupling to the field, does not allow a closed set of equations to be 
deduced. Quang and John attempted a treatment in which they assumed that no more than one photon was ever 
present in the output field [ jl4| , but in a subsequent comment this assumption was shown to be equivalent to the 
Born approximation, and the Monte Carlo wavefunction method employed to deal with this problem was shown to 
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The memory of a system with a quadratic dispersion relation in one dimension is effectively infinite. This is evident 
in the undriven case where photons can be permanently localized near the atom, which can remain partially in the 
excited state forever. Driving the atom in the presence of this dressed atom-photon bound state will cause the bound 
state to become more highly occupied. A driven atom in a PBG material is analogous to the pumped atom laser 
model, which does not have a steady state without the inclusion of a loss mechanism for the bound state p(J . In a real 
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system the build-up of energy in the bound state would reach some limit as imperfections in both the PBG material 
and the quadratic approximation became significant. In other words, describing the resonance fluorescence of an atom 
in such a material requires inclusion of the losses in order to produce a physical result. This will result in a finite 
memory time for the evolution of the atomic system. In addition, the calculations so far have dealt with an isotropic 
model of a bandgap but more appropriate for artificially created bandgap materials is the anisotropic bandgap model 
pl[ which leads to memory function that decays much more rapidly than that for an isotropic bandgap. In this paper 
we will present a method which is able to calculate the evolution of a driven atom near a PBG providing the memory 
time of the system is sufficiently short, and the coupling is sufficiently weak. 

In the next section we will describe our model for the driven atom, and demonstrate how a formal solution for the 
reduced density matrix of the atom can be generated. In Sec. Ill we will describe a numerical method for producing 
results with this formal solution. Sec. IV shows how our method reproduces the exact result from a similar, soluble 
system, and then describes the spectrum of the resonance fluorescence near a band gap in a PBG material with 
varying levels of detuning of the atomic resonance from the bandgap edge. 



II. MODEL 



We consider a two level atom with atomic energy level separation Hloq. It is described by the Pauli spin matrices tr, 
<7^ and a z , which satisfy the commutation relations [c^er] = a z , and {er^c} = 1. The atom is driven by a classical 
field of frequency ojl, and is modeled by the Hamiltonian 

H = — a z + — {exp(-iLj L t) er' + cxp(iuj L t) a}, (l) 

where fi is the strength of the driving field experienced by the atom. 

Under the rotating wave approximation [ p2[ the Hamiltonian for the atom coupled to a large number of discrete 
modes of the electromagnetic field can be written as 



H = H + [hokblbk + M9*kb{<J - g k b k ^)\ (2) 



where b k is the annihilation operator of the fc-th mode, characterized by the frequency u k , and is coupled to the atom 
with the coupling strength g k . It is important to note that the coupling is linear in the field operators. 

For the formal manipulations which we will perform below it is convenient to work in the interaction picture. In 
the interaction picture the Hamiltonian can be written as 

H I {t)=ih[g{t)*{t)-t{t)^{t)l (3) 

where the driving field [ p3| , is defined by 

m^Y,9kb k e~^\ (4) 

fc 

and the interaction picture atomic operator alt) is given by a(t) — U^it) a Ua{t). The free evolution operator Uq in 
this equation is defined by Uo(t) — e~^ Hot . 

A quantity fundamental to any attempt to describe this interaction in terms of the atomic system alone is the 
memory function, defined as the commutator of the driving field with its complex conjugate; 

/ m (M') = K(*U + (*')]■ (5) 

For t > t' , the memory function is the probability amplitude of absorbing a photon at time t that was emitted 
previously at time t' . It is therefore the Feynman propagator of a virtual photon. Non-Markovian behavior arises 
when one cannot consider the lifetime of the virtual photons to be vanishingly small. The Born-Markov approximation 
can be made for free space resonance fluorescence, as (for sufficiently high frequencies) the linear dispersion relation 
(uj k = ck) produces an infinitely narrow memory function. Unfortunately, near a photonic band gap the dispersion 
relations can be quadratic, and this leads to a very long memory time. 

The combination of non-Markovian evolution and the nonlinear behavior of a driven atom renders this problem 
intractable to either Born-Markov master equation methods applicable to free-space resonance fluorescence p9[ , or 
Laplace transform methods that have been applied to the case of spontaneous emission into a band gap material [ [To[ . 
Instead, we approach this problem by a direct numerical solution of the evolution of the reduced atomic system. 
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The first important step is the derivation of a formal equation for the reduced density matrix of the atomic system. 
The reduced density matrix is defined by tracing over the field degrees of freedom of the density matrix for the full 
interacting system. In the interaction picture it has the form 



p(t) = Trfidd \ *T 



Ptot 



T 



1 Jl<iuHi{u) 



(G) 



where (T)^T denotes (anti-)timeordering. Assuming an initial factorized state of the atom and field, ptot = P® Pficid, 
and that the field is initially in a generalized Gaussian state |24|] (e.g. a thermal state), a formal expression in terms 
of only system variables can be determined for the right-hand side of Eq. (^) . Here we will perform a derivation for 
the case of an initial vacuum field, /9fi c id = |{0})({0}|. 
We use the identity p3| 



= i T a {v + v V-} 



(7) 



where T a is a time ordering operator which denotes a time ordering on the atomic operators only, and where we have 
defined the operators 



V+ =cxp (- J ds^{s)a(s) 



Vq = cxp I — du dv f m (u,v) a\u)a(v) 



V- = exp 



ds£(s)crt(s) 



where f m {u,v) is the memory function, defined by Eq. (^|). Similarly, we can write 



T 



(8) 
(9) 
(10) 

(11) 



for the anti-timeordered term. 

In order to greatly simplify the manipulation of these time-ordered quantities we introduce a notational apparatus, 
first used by Feynman p5| ], whereby unprimed atomic operators are assumed to be time-ordered to the left of p and 
primed operators anti-timeordered to the right. We can then write Eq. (|^) as 



p(t) = Tr ficld \y+V V-{p® |{0})({0}|)^'V 'V + /f } , 



(12) 



where Vj' — Vj(£, a', cr 1 '). This notation allows us to manipulate the atomic operators as if they were c-numbers 
while keeping track of their correct ordering in the expression. We can use the cyclic properties of the trace over the 
field variables to write Eq.(|l^) in the form 



(13) 



P (t) = ({o}|y_'V 'V + 'V + %v_|{o})p. 

We now normally order the field operators in this expression by making use of the relations e 
e B e A e ^[A,B] £ or ^ w0 p era t or s A and B satisfying [A, [A, B]] = [B, [A, B]] = 0. The normally-ordered field operators 
are annihilated on the field vacuum and we are left with the result 



A+B 



e A e B e -±[A,B] 



P(t) 



e CW P, 



where 



£(<)= J du J dv {[a'\u)-^(u)]f m (u, v)a(v) + [a(u) - a'(u)] fc(u, v)a'\v)} . 



(14) 



(15) 



This is a formal solution that contains no reference to the field operators. The only reference to the particular 
characteristics of the field is via the memory function. The a^(u)a(v) and a' (u)a'\v) terms arise from the vacuum 
field to vacuum field evolution while the terms linking the two sides of the density matrix, a'^ (u)a(v) and a(u)a'\v), 
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arise from tracing over all the photons that are irretrievably emitted into the field. Note that only under the Born- 
Markov approximation [where it is possible to replace f m (t, s) by jS(t — s)] does the evolution matrix, factorize 
at each time, as in this case the double integral reduces to a single integral, ft is then possible to write down a 
ordinary differential equation for p(t). In general, however, this is not possible and it is necessary to use the formal 
solution Eq.([f|). 

The theoretical difficulty with using this formal solution is not the size of the system, which can be described by 
three real numbers, but with the implicit time ordering of the system operators, which in general leads to an infinite 
number of ordered terms. Physically, these higher order terms are the processes associated with increasing numbers 
of photons, and will become increasingly less significant in a system with a finite memory time. A driven system 
without a finite memory time may lead to unphysical behavior |L9| , so we are particularly interested in systems with 
a finite memory time. When this timescale is shorter than other dynamic timescales we are describing a Markovian 
system, and otherwise we are describing a non-Markovian system. 

In this work we consider a numerical approach to solving the evolution matrix e^K We discretize time and then 
make the approximation of a finite memory time T m . In other words, we assume that we can effectively put f m (t, s) = 
for all s < t — T m . Our algorithm considers all possible processes that can occur with these restrictions. From Eq.(|T^) 
it is clear that it is not possible to determine p(t + At) from only p(t) as at time t there are photons in the field, 
represented by the presence of the memory function that will re-interact with the atom at time t + At. These photons 
have to be taken into account in any algorithm. It turns out that under the finite memory time approximation we 
can follow the evolution of these photons by introducing a finite number of additional, "virtual", density matrices, 
p k , k 0. The evolution equation for the reduced density matrix together with the virtual density matrices has the 
form 

p k (t + At) = J2 D kj(t)P i (t) (16) 

3 

where Dkj(t) = Dkj a(t),a^(t),cr'(t),a (t) is an evolution superoperator and p°(t + At) w p(t + At), approximates 
the reduced density matrix. 

The method is based on transforming the numerical algorithm for propagating a non-Markovian differential equation 
into a Markovian form. To motivate this idea, consider the following elementary example of an integro-differential 
equation, 



dc(t) 
dt 



= -iu c(t) - / f m (t - s)c(s) 

Jt-T m 



with a finite lower-bound on the integral. This integro-differential equation can be propagated numerically via the 
algorithm 



A* 2 

z n+ i = (I - iu At)c„ — 



/oC„ + 2 ^ fn-j + lCj + fMC n -M+l 
j=n-M+2 



0{At 2 ) (17) 



where c„ = c(nAt), /„ = f m (nAt) and M + l = T m /At. We have used the trapezoidal rule for the integral and Euler's 
method to propagate the differential equation forward. Due to the finite memory time, this algorithm can be written 
in an explicitly Markovian form by introducing M — I extra variables b k . The equivalent Markovian algorithm is 

At 2 

cn+i = (1 - iojAt)c n - —f c n + 6 ? \ + 0(At 2 ) 
6i +x = b 2 n - At 2 f 2 c n 



J n+1 

b 2 l+1 = bl - At 2 f 3Cri 



h M-l At 2 
0„+l =-— jMCn 

with the initial condition £>q = 0. In this particular case, this transformation does not offer any advantage, however, 
this example is analogous to the algorithm for propagating the reduced density matrix that we will present in the 
next section. In that case, due to the presence of time-ordering, reformulation as a Markovian algorithm becomes 
invaluable. 
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III. NUMERICAL METHOD 



As a first step in determining an algorithm let us evaluate the integrals in the exponential in Eq.(|l4|) on a grid, 

j=N 

£{t)= J2 AtLj+OiAt 2 ) (18) 

3=0 



where 



L, 



]T At {(a'] a])fi- k a k + (a 3 a' r J[ k a'{) . (19) 



fc=max(j-M+l,0) 



For convenience of notation we have written tj — jAt where tjy = t, a.j — a(tj) and /„ = W n f(nAt). The weights are 
given by 

n = 

W n = { i, n = M-l , (20) 
otherwise 

as we are using the trapezoidal rule (as the time ordering produces an integrand that is effectively discontinuous at 
each time step there is good reason to believe that a higher order rule will not improve accuracy). Note that we have 
also explicitly truncated the sums after M > T m / At as the memory function is assumed to be zero for all additional 
terms. We can now write 

e m=Y[e AtL * + 0(At 2 ), 

3 

= Y[(l + AtL J )+0(At 2 ), (21) 

j 

where in the second line we have expanded out the exponentials at each time, tj, to first order in At. Note that it is 
still necessary to timeorder the operators in this expression. Let us attempt to write this expression in terms of the 
above mentioned density matrices p k at one particular time n — 1, where we choose N — M > n > M as this is the 
most common case. First split the evolution matrix at n + M and n, 

N n+M-1 

e m = ]J (...) l[ (l + AtL^p + OiAt 2 ), (22) 

k— n+M j—n 

where it is convenient to separate out the factor p = Y\"Zq(1 + AtLj)p(0) as it does not take part in our subsequent 
manipulations. Then split the sums over k in Lj for j < n + M into parts containing interaction picture operators 
before and after n, 

(1 + AtLj ) = Xj , n + At(a'l - a])y jtn + At(aj - trj>i,n, (23) 

where 

j 

x jtn - 1 + At ]T At [{a'] - a])f^ k a k + (a, - ^)f*_ k a'l\ (24) 

k—n 

n-1 

Vi<n = ^2 ^tfj-kCTk (25) 

k=j-M+l 
n-1 

zj, n = Yl A */;_ fc <4. (26) 

k=j—M+l 

Consider the M factors from j = n to j = n + M — 1. Due to the truncation of the sums at M terms, it is only these 
factors that contain interaction picture operators at time n. We can expand this product of M factors, each of which 
has three terms, as a sum over 3 M terms 
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n+M-1 3 JU ~-1 

(27) 

j—n m—0 

where g m ,n contains only operators after n — 1 and h m>n contains only operators before n. Note that since y n +M-i,n — 
z n +M-i,n = we only need to consider a sum over 3 M_1 terms. To write down the explicit form for g m ,n and h m ,n 
it is convenient to introduce a new notation to label the terms. We write the index m in trinary (i.e. base three), 
where the fc th trit (bit in base three) tells us which of the three factors to choose from the (1 + AtLk) term (shown 
in Eq. (p3|). Our new index then lists the trits which correspond to the first term, {Ai, ■ ■ ■ ,A p }; the second term, 
{Bi, • ■ • , B q }; or the third term, {C\, • • • , Cu-q-p}- These lists replace the m index. In this notation we can write 

P <? 

9n — x n+M-l.n [ J_ X n+Aj _!,„ j| A.t(CT n+Bk -l ~ a n+B k -l) 

j=l k=l 

M—q—p 

x J] At (°n +Cl -i - <r' n+Cl -i) (28) 
1=1 

and 

q M—q—p 

h^ Bl ''"' B "' Cl '"'' CM -"-^ = y n +B k -l,n it Z n +Cl -l,n. (29) 

fe=l Z=l 

where, if the lists {Bi, • • • , B q ; C\, • • • , C^_g_ p } are the empty set, 0, by construction we have = 1. This allows 
us to define a set of 3 M_1 density matrices at time n — 1 of the form 

)) {b 1 ,..,b,;C 1i -.,c m -,_ p } = /,;/-• ,,/;,:-• ,^.r, ; . ,- /K (30) 

Note that p\_ x is defined such that it corresponds to p(t n -\) up to order At 2 . 

If we can determine the transformation from p„_x to p 3 n we will have achieved our aim of writing the non-Markovian 
evolution in terms of a Markovian evolution in an extended space. This is done by writing p> n in terms of p^-i- Noting 
that yj,n+i = Vj,n + Atfj-ndn and Zj, n +i = z 3 , n + At}*_ n a'\ we can write 

q M — q—p 

p {Bu...B t ;0 lt ...Cu- t - p } = ~[[lyn + B k ,n + At/ Bfc (7 n ] j\ [ Z n+C 3 ,i + Atf c .a'l] 

k=l j=l 

x [x n>n + At(a' n - cr^)y ni „ + At (er„ - a' n )z nin ]p. (31) 

Expanding out the first two sets of factors and keeping only first-order terms the right hand side can now be written 
in terms of the density matrices at n — 1, 



n {Bi,"-,Bq;Ci,---,CM—q—p} _ 
Pn 



l + At (a n - a-Dfoan + At (<7„ - cJ/ V n 
q 

{B 1 +l,--;B k _ 1 +l,B k+ - L +l,~-,B q +l;C 1 +l,--;C M -<i-p+l} 



{Bi+1,-,B,+1;C i +1,-,Cm-,- p +1} 
Pn-1 



+ AtC7 n J2f Bk p n Ji 



k=l 
M-q-p 

, A , it \ - f * {Si + l,--- : S, + l;Ci + l,---,C fe _i + l,C fc + i + l,---,CM-<j-p + l} 
"T^" 7 n 2-1 JC k Pn-l 
k=l 

- /T \\ n {1-B 1 + l,-,B q + l;C 1 + l,-.,CM- q - P + l} 

v n a n)Pn-l 

+At(a n - ^) p {^+i.-,B,+i;i,c 1 +i,-,<7 M _ 4 _,+i} + 0(A ^ (32) 

Physically, during a short enough interval only one photon can be emitted or absorbed at a time. The higher order 
terms in At correspond to two or more particle emissions or absorptions during this interval and hence we have 
dropped them. That is, except for the terms At 2 (a' r ^ — a} l )focr n + At 2 (a n — a' n ) /q a'^ , which we include because often 
the memory function diverges at the origin and /o ~ 1/Ai, so this term is effectively of order At. The additional 
density matrices are virtual in the sense that they represent states that are never observed in the system behavior 
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and also in the related sense that they follow the paths of virtual photons. This equation has the form of Eq.(|lJ) and 
therefore the right hand side can be put in the form of a superoperator acting on a density matrix. In this way it is 
possible to propagate the non-Markovian evolution forward in time. 
Any expectation values are determined from p^ n alone; 

(a n ) = Tr{a„^}. (33) 

For this reason it is convenient to introduce the projection operator P defined by p\ = Pp , where p\ is projected 
out of the "vector" of matrices p n which includes the reduced density matrix and all of the virtual density matrices. 
In this notation we can write the expectation value Eq.(p3) as 



(a n ) = Tr{a n PpJ (34) 

and the evolution equation (5^) as 

P n = D nP n _ V (35) 

where D n = D(a n , <jt , a' n , a' n ). Often one is not only interested in expectation values at a single time but also in 
temporal correlations. For example, the correlation function (with N2 > N±) 

(^ N2 (J Nl ) = Trjerj^gjVs.JVi} (36) 

where, assuming that the initial state of the field is the vacuum state, we have written 

QN 2 , Nl =<T Nl e^ N ^p. (37) 

In general, techniques for calculating the correlation function such as the quantum regression theorem (see, for 
example, j^6|) will not work in the non-Markovian case as there is no factorization of the evolution matrix e ci * N2 ^ 
at time N%At. However, it is not difficult to see that this type of correlation function is easily calculated numerically 
within the above formalism. One simply propagates the set of density matrices p forward from n = as usual; then 
at n = N\ multiply by <Jn 1 ', taking this as the new initial state one can propagate up to time N2 giving qn 2 ,n 1 to the 
accuracy of the algorithm. In symbolic notation this procedure is written as 



-JVa < — Ni 



Qn 2 ,N! 



p n fe ^ 1+1 ^ru^o> <») 



where denotes a time-ordered product. 

For numerical simulations it is often more useful to work in the Schrodinger picture. Equation (|3^) can easily 
be put in the Schrodinger picture because the operators are all acting at the same time. We simply multiply by 
Uo(nAt)UQ (nAt) which transforms the density matrices at time n to the Schrodinger picture and then write the 
interaction picture operators on the right hand side in terms of the Schrodinger picture operators. The resulting 
equation has the same form as Eq.(|32|) (with the density matrices in the Schrodinger picture) except the right hand 
side is multiplied by U (At)U^ (At); 

£ = D(a,*\a'^)U Q (At)U'J{At)£_ x . (39) 

The advantage of the Schrodinger picture is that the superoperator D does not change with time and so only needs 
to be constructed once at the start of the simulation. 

It is worth noting that the above superoperator representation does not need to be abandoned in numerical simu- 
lations. We can, for example, flatten the density matrix into a column vector so that both the primed and unprimed 
operators (which consequently become 4x4 matrices rather than 2x2 matrices) act only on the left-hand side of 
this flattened density matrix [^tJ. The right-hand side of Eq. (|39|) then represents the 4(3 M ) x 4(3 M ) sparse matrix 
D(a, a', a\ a'^) U (At)U' Q \At) multiplying the (3 M )4 x 1 column vector p^. Note that the initial vector p Q has all 

elements zero except for = p. Also note that we have written M not M — 1, as some virtual density matrices are 
initialized at each time step. 
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IV. RESULTS 



Before we present the results of numerical simulations using the above algorithm, it is important to make some 
general comments about the limitations of the algorithm. As mentioned above, to propagate the non-Markovian 
evolution requires a vector of size (3 M )4 x 1 to be updated at each time step, thus computer memory considerations 
put a severe limitation on the number of time steps per memory time, M, that the algorithm can deal with. In our 
case, we created the time independent evolution matrix at the beginning and were in fact limited by the number of 
non-zero entries of this matrix that could be stored in computer memory. On a computer with 128MB of RAM we 
were limited to a maximum of M = 11. This obviously poses quite a problem when one wishes to consider quite 
non-Markovian situations but also require high accuracy. However, for test cases where the system was not strongly 
non-Markovian, the results of our simulations show good agreement with calculations using standard methods of 
analysis. 

In the rest of this article we consider the case of on- resonance driving (loq — lul)- In a frame rotating at the atomic 
frequency the Hamiltonian for the atom reduces to 

Ml , i . 
H Q = — + a) . 

This Hamiltonian gives rise to the following behavior 

o-(t) = eTt Hot <je-K Hat 

= l(<r x + \+){-\e mt -\-){+\e-^) 

where a x = cr + er and the states |±) are defined by cr x |±) = ±|±). When the atom is coupled to a propagating 
field these three terms give rise to three distinct peaks in the spectrum at u> — luq — 0, ±f2. This is referred to as the 
resonance fluorescence triplet [ p8|]29| ] . 

With a structured electromagnetic field, these emission peaks may be altered. In the following sections we use our 
algorithm to consider the case when the atom is in a cavity tuned to one of the side peaks and the case when one of 
the side-peaks is at a frequency where no photons propagate in a bandgap material. 



A. Atom in a Cavity 



In order to demonstrate the validity of our algorithm we consider a simple system that can also be described by 
a Born-Markov master equation (and so can be solved by traditional methods) and compare the results of the two 
methods. The system that we consider is an atom coupled to a single mode of a cavity which is in turn coupled 
to a continuum of propagating modes of the electromagnetic field via a partially transmissive mirror. The standard 
method of finding the evolution of this system is to treat the atom and the cavity as a system and the external field as 
a reservoir which can be traced over. This produces a Born-Markov master equation for the atom and cavity modes. 
This can then be solved by standard methods p6| . One can also make a different choice of system and reservoir; 
consider the evolution of the atom alone while treating both the cavity and the external field as the reservoir; giving 
rise to a non-Markovian problem. 

The usual description of this system is in terms of the Hamiltonian 
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Hx = H + H + ih^J [aV - a^a] , (40) 
where 

H = Twa)a + h[ duj I utf (<j)b(u) + i-£= [tf{w)a - 6(cj)a f ] ) , (41) 

J-OO I V27T J 

and where a is the annihilation operator of the cavity mode, b(u>) is the annihilation operator for the external 
propagating electromagnetic field mode of energy Tko, a is the lowering operator for the atom, and Hq acts only on 
the atomic operators. 

Choosing the atom as the system and eliminating the cavity mode it is possible to describe this system by the 
alternative Hamiltonian (see Appendix. ^ ) 
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H 2 =H Q + h oho {ujJ(oj)c{uj)+i[\*{uj)J(oj)a-\(uj)c(uj)^}}, (42) 

J —OO 

where c(u) is a new bose annihilation operator and the, now frequency dependent, coupling constant is 

Ko>) = \[3zt, K —^- (43) 



2tt »( w - v) - f 
The memory function in this case is therefore given by 



oo poo 



— oo •/ — oo 



fm(t-s)= duj I duj 1 A*(w')A(w)[c(w),cT(w')]e" Mt+ws (44) 



X / ~ duj K l e- 1 ^-^ (45) 



2tt 



H 2 



= ie xp^-iu(t-s)-—\t-s\j. (46) 

Transforming to a frame rotating at the atomic frequency appends a factor e luJ o{t-s) ^ Q mem0 ry function. This 
model differs from that of the bandgap only by the functional form of the memory function. 

We compare our results for the system described by H 2 using the non-Markovian algorithm introduced in Sec. 
H| with the solution of the Born-Markov master equation for the atom plus cavity system as described by Hi . In 
Fig.|l| we plot the time evolution of the probability for the atom to be in the excited state. The internal spectrum, 
S((JT, defined as the Fourier transform of the steady-state correlation function (<T^(r)cr) ss — (o^)ss(c)ss , is plotted in 
Fig.H showing a significant asymmetry in the height of the side peaks due to the tuning of the cavity to only one 
side band of the resonance fluorescence triplet. This is a completely non-Markovian effect as it arrises directly from 
the structured nature of the field, i.e., the fact that some frequencies of the field are easier to emit into than others. 
Both the probability and the spectrum show good agreement with the results of the traditional method (we used the 
numerical routines in the quantum optics toolbox |27[]). 



B. A Bandgap Material 

In this section we apply our algorithm to resonance fluorescence in a band gap material; a problem that cannot be 
solved by any traditional methods. Following the discussion of Vats and John |2lf| we consider an anisotropic bandgap 
model under the effective mass approximation. This approximation is appropriate to fabricated band gap materials. 
In this model the band edge is associated with a specific point in fc-space, k = ko, and the dispersion relation close 
to the upper edge of the gap has the form 

?i(k-k ) 2 , s 

= u g + K 2 ' ■ (47) 

where /i is called the effective mass in analogy with the dispersion relations of massive particles. 

This dispersion relation gives rise to a memory function with the asymptotic behavior f m (r) ~ r~ 3 / 2 so that, in 
contrast to an isotropic bandgap (where fra{r) ~ r -1 / 2 ) we are able to define a finite memory time. 

The memory function given by / m ( r ) °^ r~ 3 / 2 suffers from a singularity at r = 0. This is due to the absence of a 
high-frequency cutoff, which must be included to model a physical system. We include this high frequency cutoff by 
considering the slightly modified memory function, 



_i[5T+ir/4] 

(1 + Ar) 



where (shifting to a frame rotating at the atomic frequency) 6 = luq — uj g is the detuning between the band edge and the 
atomic resonance frequency and is a measure of the coupling strength between the atom and the electromagnetic 
field. The parameter A is infinite when there is no upper frequency limit, so for large A this memory function is 
effectively the same as that derived by Vats and John. It is useful to introduce the quantity 







/•OO 

7 = 23?/ dr / m (r)| 5=0 
Jo 



which corresponds to the Born-Markov damping rate for this memory function. We use this quantity as a scale against 
which we can measure the various rates. 

For large A, the factor (1 + At) 3 ^ 2 has a sharply rising behavior as r — > on a time scale much faster than the 
system variables. This is a common feature of physical memory functions and requires special treatment in numerical 
calculations. We generalize the trapezoidal rule used for the atom in a cavity case by calculating weights that depend 
on the first two moments of (1 + At) 3//2 and derive an extended two-point integration scheme on a uniform mesh as 
described in Ref. |M] . Explicitly, we used the weights 

Wq — w®, n = 

W n = { — (M - 2)wf- 2 + wf ~ 2 , ri = M - 1 

(n + 1)wq - - (n- 1)wq 1 + otherwise 



where 



Since 



J AP 



(n+l)At 

dr r j (1 + At)~ 3/2 . 



A/ 



n=M-l r (M-l)At 

dr (l + ATp 3/2 . 



E ^ = / 



this integration scheme treats the semi-singular part of the memory function exactly and allows us to treat a range 
of non-Markovian situations, including the Markovian limit, within the same algorithm. As mentioned previously we 
do not consider schemes of higher order than the two-point scheme as the numerical difficulty of evolving the system 
a single timestep increases dramatically due to the combinatorics introduced by the time-ordering. 

Let us first consider the undriven case where an initially excited atom simply undergoes decay inside a bandgap 
material. Since this case can be treated by an alternative method it serves as another test of the algorithm. In the 
case of simple decay there is at most only one photon in the electomagnetic field so we can write the state of the atom 
and field as 



|*(t))=a(i)|l)|{0})+^c fc (t)4|0)|{0}), 



where |1) and |0) are the excited and ground states of the atom and |{0}) is the field vacuum. It is not difficult to 
establish that 



da(t) _ r ' 
dt ~ 



ds f m (t - s)a(s) 



n 



with f m (t — s) given by Eq.(|4^). This integro-differential equation can be numerically simulated by, for example, 
the routine (|17]). In actual fact we used the weighting scheme described above to approximate the integral. The 
probability of the atom being in the excited state is given by P(t) = \a(t)\ 2 . In Fig. || we have compared the results 
of our algorithm to a direct numerical integration of the above integro-differential equation; the two methods show 
good agreement. In contrast to the isotropic band gap case, previously treated analytically in Ref. [0, the decay is 
always exponential. This can be seen more clearly in the corresponding plot of the logarithm of P(t), Fig. [|. Note 
that in this plot the discrepancy between the two methods is apparent. Our algorithm gets more accurate as A is 
increased as this allows a smaller stepsize to be taken. These plots were made with the parameter value A = 3 x 10 2 7 
and the variation of the damping rate with detuning from the bandgap is evident. However, when this parameter is 
increased the situation becomes more and more Markovian, despite the long time tail of the memory function. For 
A = 10 5 7 there is no noticeable difference from the Markovian case; the decay rate is exactly that predicted by the 
Born-Markov result, and the detuning becomes irrelevant. This is shown in Fig.||. Note that the S < case requires 
the driving laser to be propagating in a direction that is perpendicular to the direction of ko for the anisotropic 
bandgap. 
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For non-zero driving there are no traditional methods available and the results can only be determined by a quantum 
non-Markovian treatment. In Fig. [| we have plotted the probability amplitude for the atom to be in its excited state. 
The driving induces oscillations in this probability and the atom tends to a steady state close to equal probability of 
being in the excited state and in the ground state due to the reasonably large driving. This behavior is qualitatively 
equivalent to the Markovian case. The steady-state correlation function, C(r) = (a^(r)a) ss — (<t^) ss (ct) ss , is plotted 
in Fig.Q. Taking the Fourier transform of the steady-state correlation function, C(r), gives the so called internal 
spectrum of the atom, S(oj). As in the atom in a cavity case this quality allows us to clearly see the effects of the 
non-trivial structure in the field on the behavior of the atom. In Fig.|^ we have plotted the internal spectrum with 
A = 3 x 10 2 7 for various detunings from the band edge. A definite asymmetry arises in the magnitude of the side 
peaks directly due to the fact that emission at frequencies inside the gap is prohibited. Note that as A is increased this 
effect becomes less and less until the spectrum becomes equivalent to the Markovian case. In Fig we have plotted 
the internal spectrum with the same detunings but with A = 10 5 7; there is no evidence of asymmetry and the values 
of detuning become irrelevant just as in the case of the undriven atom. 



V. CONCLUSIONS 



A number of different methods of dealing with non-Markovian systems have appeared in the literature recently. 
These methods can be broadly separated into two main categories: methods that involve extending the system size 
by the addition of real or imaginary modes p^-p3] such that the new extended system evolves via a Born-Markov 



master equation, and methods that attempt to evolve the non-Markovian evolution directly [|23| , |34 - 37 1 ■ Note that the 
approach in Ref. ||(| takes advantage of the special case when the free evolution of the system has an analytical solution 
to derive a Markov equation for the system. Conceptually, the present work falls into the last category, however, in 
order to propagate the non-Markovian evolution we have effectively extended the system size by introducing virtual 
density matrices that follow the paths of the photons in the field before they are either irretrievably emitted or are 
reabsorbed by the atom. The two approaches have very different origins and the way in which the extended system 
is constructed is fundamentally different. The present method seems to be slightly more general as no restrictions are 
made on the form of the memory function, however, time discretization plays a much more fundamental role in our 
algorithm so it would be very interesting to see on what level these two approaches are related. 

The strength of the present algorithm lies in its ability to accurately model weakly non-Markovian systems. An 
appropriate application would be to determining the validity of the Born-Markov approximation in borderline cases. 
The biggest constraint on the algorithm is the power law increase in the number of virtual density matrices that 
need to be stored to propagate the non-Markovian evolution. This is a fundamental constraint and arises from 
the inherent difficulty in numerically simulating quantum mechanical systems in general. Note that the algorithm 
presented here was derived without using the algebraic relations of the atomic operators. It is therefore applicable to 
other system-bath couplings that are linear in the bath operators. 

With our algorithm we have determined the evolution of a driven atom radiating into a bandgap material with an 
anisotropic bandgap; a situation that cannot be treated by any traditional methods. The asymmetry in the internal 
spectrum of the atom clearly demonstrates the effect of the non-trivial structure of the field on the behavior of the 
atom. This behavior cannot occur under the Born-Markov approximation. In addition, we have demonstrated the 
transition to the Born-Markov limit, which occurs despite the presence of a long time tail on the memory function. 

In this article we have presented a density matrix version of the algorithm. In this case the system size was so 
small that a density matrix treatment was appropriate, however, a wave function version can also be derived in a 
similar fashion. The wavefunction version is applicable to simulating non-Markovian quantum trajectories that also 



solve for the reduced density matrix |E3j,p4 35 37]]. In the non-Markovian case increasing system size is obviously an 



even more difficult constraint to numerical commutation than in the Markovian case so wave function methods will 
be important. 
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APPENDIX A: NON-MARKOVIAN MODEL FOR AN ATOM IN A CAVITY 



This appendix shows how the system in section IV A can be written in terms of a non-Markovian model. 

Equation (^) is valid in the high-frequency limit |36|]: it assumes the rotating wave approximation and that the 
fundamental frequency v is much greater than the spread of frequencies about the fundamental frequency such that 
we can replace the lower bound of the integral over the bath modes by — oo. Note that as we are considering a one 
dimensional field in the continuum limit it is convenient to work in terms of the frequencies rather than the mode 
labels. 

The proof that the two Hamiltonians (f40| ) and ( ff2| ) are equivalent consists in demonstrating that they are related 
by a unitary transformation. Let the initial time of interest be t — 0. For the moment let us restrict our consideration 
to the cavity mode-field system described by the Hamiltonian H [see Eq.([|l])]. Consider the following scenario: let 
the coupling between the cavity and the field be time dependent and let it increase slowly such that at a long time in 
the past (t < 0) the coupling is turned off and at the initial time (t = 0) it is equal to k, 

*>S- (A1) 

where e is a small positive number. The adiabatic limit of an infinitely slow switching on of the coupling is given by 
the limit where e € * — > as t — > — oo. and e — > The Heisenberg equations of motion for the cavity and field modes 
for t < are given by 

^ =-i [2(t), H(t)] = (-i v - e )a(t) - -£= f dw b(u,t) (A2) 



dt Ti l V2^ J- 

= -i[6( w , t),S(t)] = -iub{u) + -|=2(t), (A3) 

where we have simplified the equations by making the replacement a(t) = a(t)e et . Note that if limt ^_oo a(to) is finite 
this implies that 

lim a(t ) = 0. (A4) 

to — > — oo 



Solving Eq.(A3) formally for b(u),t) in terms of b(uj,t ), where to < t, and substituting this solution into Eq. (A2) we 
derive 



da(t) 
dt ~ ' 

where we have used the relation 



(2 \ roo 
iv + e + ^-)a{t)--j= / duj b{uj,t )e' l ^ to \ (A5) 
2 / V27T J-oo 



1 f°° 

5(t-s) = — duj e -^ ( *- s) . (A6) 

27r J-oo 

Solving this formally for a(0) in terms of a(<o) we find 

!, K I" 00 ~ 1 - „[i(^-w)+e+^-]t 

5(0) = d(t )e^ + ^ j= J du; b(u;,to)jr, ~ ^t - (A7) 

V27TJ-0O [i{v - uj) + e + %J 

where b(u),to) = &(w,to)e Iwto - Letting to — > — oo and e — > (and assuming that b(u),to) exists in this limit) this 
solution becomes 

o(0) = 2(0) = -= / dw b(Lj,-oo)— — (A8) 

As the system contains no bound states the Hamiltonian in the same limit is given by 

/oo 
duj wtf(tJ,-oo)b(u,-oo). (A9) 
-oo 

Therefore, the unitary transformation given by 



12 



U = lim 

tQ— 00,€— *0 



cxp < — / ds H(s) 



(A10) 



diagonalizes the Hamiltonian H(0) and gives a(0) in terms of the diagonalizing modes. From Eq.(A4), (A.8) and (AS) 
we see that this is just the transform required to write the Hamiltonian Eq. (f40|) as Eq.(|42|) with c(uj) — b(oj, — oo), 
thus demonstrating their physical equivalence. 
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FIG. 1. Time evolution of P(t) from an initial excitted state for the parameters: Q = 47, coo — v = Q, and k 2 = 87. For the 
non-Markovian algorithm we used a memory time of M — 11 and a time step of At — 1/14. The solid line is the result of the 
non-Markovian algorithm and the crosses the result of a numerical solution of the master equation for the extended system of 
atom plus cavity mode. 
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FIG. 2. Internal spectrum, S(ui), for the same parameters as Fig.l. Again the solid line is the result of the non-Markovian 
algorithm and the crosses the results of a numerical solution of the master equation for the extended system. 

FIG. 3. Plot of P(t) for an undriven atom initially in an excitted state for different detunings from the band edge. Parameter 
values are A = 3 x 10 2 7, M = 11 and At = 1/50. The various detunings in the plot are 8 = IO7 (solid line), 8 — (dashed 
line) and 8 = — IO7 (dotted line). 



FIG. 4. Same as Fig. 3 but with P(t) on a logrithmic scale. The straight lines indicate no deviation from exponential decay. 

FIG. 5. Plot of P(t) for A = 10 5 7 with other parameters the same as Fig. 3. The lines corresponding to different detunings 
are indistinguishable. 

FIG. 6. Plot of P(t) for a driven atom for the parameter values: Q = IO7, A = 3 x 10 2 7, M = 11, At = 1/50 and 8 = IO7. 



FIG. 7. Steady state correlation function, C(r) = (a'(T)a) ss — (a') ss (a) ss , of the atom with the same parameters as Fig.6. 

FIG. 8. Internal spectrum, S(u), for £1 = IO7, A = 3 x 10 2 7, M = 11, and At = 1/50. The various spectra at different 
detunings have been displaced for ease of comparison. The detunings are 8 = IO7 (solid line), 8 = (dashed line) and 8 = — IO7 
(dotted line). 

FIG. 9. Plot of the internal spectrum for A = 10 5 7 with other parameters the same as Fig. 8. The lines corresponding to 
different detunings are indistinguishable. 
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